function  results =  bond_pricing_solve_boundary_of_w( w,  par,  model_sol,  options )
% This function solves v(w, lambda) as a function lambda for either w=0, or w=1. 
% v_fitted_lambda is the function v along the lambda dimension, for the fixed w.
% Other values
if(w~=0&&w~=1)
    disp( 'Error in bond_pricing_solve_boundary_of_w: w must be either 0 or 1!' )
end

updating_dampenning = 0.8;

pi=par.pi;   tau=par.tau;
lambda_star = par.lambda_star;  hat_kappa0=par.hat_kappa0;
if( isfield(options, 'safe_bonds') && options.safe_bonds ) % Safe bonds 
    hat_kappa0 = 0;   par.hat_kappa0=0;   % Note: it is important to change par_bonds.hat_kappa0, since we will call the function "bond_pricing_ode_lambda", which needs this parameter.  However, note that this will only change the values of this struct in the internal environment of this function.
end
if(w==0)
    rf = model_sol.rf0;
elseif(w==1)
    rf = model_sol.rf1;
end
kappab = model_sol.kappab_fitted_matrix(w,lambda_star);
lambda_star_post_jump = lambda_star + par.kappa_lambda_fun(lambda_star);
lambda_tau =  1/tau - lambda_star*pi; 
v_lambda_star = ( lambda_star*pi/(1-kappab)*(1-hat_kappa0)+lambda_tau )/(lambda_star*pi/(1-kappab)+ rf +lambda_tau);
v_new = @(x) v_lambda_star*ones(size(x));    

if( par.model_name == "benchmark"  )
    diff_vec = 0;
    lambda_vec = par.lambda_grid;
else
    diff = 1;  diff_vec=[]; iter_max=100;  iter=0;
    while( diff>0.0001 && iter<iter_max )
        iter = iter+1; 
    %     dvdlambda = @(lambda,v) ( rf0*v + lambda + lambda*pi*hat_kappa0 )/ par.mu_lambda_fun(lambda);
        v_last  = v_new;
        v_lambda_star = ( lambda_star*pi/(1-kappab)*(1-hat_kappa0)+lambda_tau - lambda_star*(1-pi)/(1-kappab)*( v_last(lambda_star) - v_last(lambda_star_post_jump) ) )/(lambda_star*pi/(1-kappab)+ rf +lambda_tau);
        dvdlambda =  @(lambda, v) bond_pricing_ode_lambda(  lambda, v, par,  model_sol,  w,  v_last );
        [ lambda_next, v_next ] = ode45(  dvdlambda, [lambda_star+1e-10, par.lambda_H ],  v_lambda_star );
%         options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot);
        [ lambda_before, v_before ] = ode45(  dvdlambda, [lambda_star-1e-10, par.lambda_L ],  v_lambda_star);
        lambda_vec = [flip(lambda_before); lambda_next(2:end) ];  % Note: this lambda_vec may change across different iteration rounds.
        v_vec = [flip(v_before); v_next(2:end) ]; 
        index = lambda_vec<lambda_star-1e-10 |  lambda_vec>lambda_star+1e-10;
        fitted_fun =  griddedInterpolant(lambda_vec(index) , v_vec(index));
        v_vec=  fitted_fun(lambda_vec);
        v_new = griddedInterpolant( lambda_vec,  (1-updating_dampenning)*v_vec +  updating_dampenning*v_last(lambda_vec)  );
        diff = sum(  abs(  v_new(par.lambda_grid) - v_last(par.lambda_grid)  )  );
        diff_vec = [ diff_vec,  diff ];
        if( isfield(options, 'plot') && options.plot) 
            plot( lambda_vec, v_vec  ); pause(2); close all;
            pause(1);   close all;
        end
    end
end

results.diff_vec = diff_vec;
results.lambda_vec = lambda_vec;
results.v_fun_lambda = v_new;
results.v_vec = v_new(lambda_vec);












